嘿嘿嘿
圓規正轉:
一樣AI slob
Day 24 把 force-directed layout 拆成兩種主要工作：
repulsion
每個 node 對其他所有 nodes
work = N * N
attraction
每條 edge 連接兩個 endpoints
work = E
Naive GPU 版本可以讓一個 thread 負責一個 node：
int i = blockIdx.x * blockDim.x + threadIdx.x;
float px = pos_x[i];
float py = pos_y[i];
float fx = 0.0f;
float fy = 0.0f;
for (int j = 0; j < num_nodes; ++j) {
float dx = px - pos_x[j];
float dy = py - pos_y[j];
float dist_sq = dx * dx + dy * dy;
if (dist_sq > 1e-8f) {
float inv = 1.0f / dist_sq;
fx += inv * dx;
fy += inv * dy;
}
}
這已經有平行度，卻留下了一筆很大的記憶體帳。
假設一個 block 有 256 個 threads。它們各自負責不同的 owner node，卻會依序讀取完全相同的 partner positions：
thread 0: pos[0], pos[1], pos[2], ...
thread 1: pos[0], pos[1], pos[2], ...
thread 2: pos[0], pos[1], pos[2], ...
...
thread 255: pos[0], pos[1], pos[2], ...
今天要解的問題不是把 N * N 假裝變成比較低的複雜度。Exact repulsion 的每一對 nodes 都仍然要算。
我們要做的是：
同一批 partner positions
只從 global memory 合作載入一次，
再讓整個 block 重複使用。
這就是 shared-memory tiling。
Day 24 用 Fruchterman-Reingold，也就是 FR，建立容易理解、可以直接執行的 baseline。原專案的 production layout 則是對齊 OGDF FMMM 的 force model。
兩者都有 attraction、repulsion 與逐輪更新，但公式不完全一樣。今天分析的 CUDA kernel 使用：
repulsion from j to i:
dx = x[i] - x[j]
dy = y[i] - y[j]
dist_sq = dx * dx + dy * dy
F_rep_x[i] += dx / dist_sq
F_rep_y[i] += dy / dist_sq
edge (u, v) 的 attraction 使用 edge 的理想長度 L：
dx = x[v] - x[u]
dy = y[v] - y[u]
d = sqrt(dx * dx + dy * dy)
scalar = log2(d / L) * d / (L * L * L)
F_attr_x[u] += scalar * dx
F_attr_y[u] += scalar * dy
最後把兩種力合併，並限制單輪移動距離：
force = (spring_strength * attraction
+ repulsion_strength * repulsion) * avg_L * avg_L
norm = length(force)
step = min(norm * force_scaling_factor, max_radius) / norm
displacement = step * force
這裡刻意保留公式差異。GPU 最佳化可以搬用相同的 tiling、memory placement 與 kernel fusion 思路，但不能在換了 force model 後，仍宣稱和 Day 24 的 FR 是逐值相同的演算法。
一個 position 有 x、y 兩個 float：
bytes_per_position = 2 * 4 = 8 bytes
現在取 256 個 partner nodes 當成一個 tile，並讓一個 256-thread block 計算 256 個 owners。
Naive 寫法中，每個 owner 都要請求整個 tile：
naive requested bytes
= 256 owners * 256 partners * 8 bytes
= 524,288 bytes
= 512 KiB
Tiled 寫法讓 256 個 threads 合作載入：
tiled requested global bytes
= 256 partners * 8 bytes
= 2,048 bytes
= 2 KiB
從程式發出的 global-load 請求來看，同一個 block 對這批 partner positions 的重複讀取最多可以從 256 份降成 1 份：
request reduction = 512 KiB / 2 KiB = 256x
但這不是 kernel 的 256 倍 speedup。
原因包括：
GPU cache 可能已經吸收一部分重複讀取
每一對 nodes 的算術仍然存在
shared memory 也有 load 指令
每個 tile 需要 synchronization
register 與 shared-memory 用量會影響 occupancy
attraction、更新座標與 launch overhead 沒有一起縮小 256 倍
所以這個 256x 只能叫做「特定 tile 內的理論 global-load request reduction」。它不能被寫成 benchmark 結果。
原專案把 tile size 與 block size 都設成 256：
#define FMMM_TILE_K 256
#define FMMM_TILE_BLOCK 256
一輪 repulsion 的資料流如下：
global pos_x / pos_y
|
| 256 threads cooperative load
v
+-----------------------------+
| shared s_tile[0 ... 255] |
+-----------------------------+
| | |
v v v
owner 0 owner 1 owner 2 ...
register register register
| | |
+--- accumulate forces ---+
核心結構是：
__shared__ float2 s_tile[FMMM_TILE_K];
for (int base = 0; base < num_nodes; base += FMMM_TILE_K) {
int nload = min(FMMM_TILE_K, num_nodes - base);
for (int k = threadIdx.x; k < nload; k += blockDim.x) {
s_tile[k] = make_float2(
pos_x[base + k],
pos_y[base + k]
);
}
__syncthreads();
accumulate_my_owners(s_tile, nload);
__syncthreads();
}
第一個 __syncthreads() 的意思是：
所有人都放完資料後，任何人才能開始讀 tile
第二個 __syncthreads() 的意思是：
所有人都用完目前的 tile 後，才可以覆寫它
少了第一個 barrier，thread 可能讀到尚未載入的欄位。少了第二個 barrier，速度較快的 thread 可能先把下一個 tile 寫進來，破壞其他 thread 還在使用的資料。
這兩個 barrier 都是 correctness 的一部分。
float2？Global memory 裡的座標分成兩條陣列：
pos_x = [x0, x1, x2, x3, ...]
pos_y = [y0, y1, y2, y3, ...]
這叫 Structure of Arrays，縮寫 SoA。
當相鄰 threads 載入相鄰索引時：
thread 0 -> pos_x[base + 0]
thread 1 -> pos_x[base + 1]
thread 2 -> pos_x[base + 2]
...
這些位址是連續的，GPU 可以把它們合併成較少的 memory transactions，也就是 coalesced access。pos_y 也是一樣。
載入 shared memory 後，kernel 把同一個 node 的 x、y 組成 float2：
s_tile[k] = { pos_x[base + k], pos_y[base + k] }
repulsion 每次一定同時使用 partner 的 x 與 y。shared memory 使用 float2，讓程式更自然地取出一個完整 position。
這是一個混合配置：
global memory
SoA: float pos_x[], float pos_y[]
目的：連續、容易 coalesce
shared memory
AoS-like: float2 s_tile[]
目的：一起消費同一個 node 的 x/y
SoA 或 AoS 沒有永遠正確的答案。要看資料在那一層由哪些 threads 載入，以及計算時會一起使用哪些欄位。
每個 owner 的資料不需要放進 shared memory。
在 tile loop 期間，owner 的座標與累加結果會被反覆使用：
owner position: qx, qy
repulsion accumulator: frx, fry
previous displacement: lx, ly
它們最適合留在 thread-private registers：
float qx[FMMM_TILE_OWN];
float qy[FMMM_TILE_OWN];
float frx[FMMM_TILE_OWN];
float fry[FMMM_TILE_OWN];
float lx[FMMM_TILE_OWN];
float ly[FMMM_TILE_OWN];
完整的資料放置可以畫成：
registers, private to one thread
owner indices
owner positions
accumulated repulsion
previous displacement
shared memory, shared by one block
current remote position tile
global memory, visible to the grid
all positions
CSR adjacency
displacement arrays
iteration outputs
判斷資料該放哪裡時，可以問：
只有一個 thread 高頻使用？
-> register
同一個 block 的很多 threads 會重用？
-> shared memory
資料太大，或需要跨 blocks 共享？
-> global memory
把所有資料都塞進 shared memory 並不會自動變快。CSR adjacency 大小隨 E 增長，而且每個 owner 只讀自己的鄰居範圍；把整份 CSR 搬進 shared memory，容量不實際，重用率也不高。
Production tile kernel 設定：
#define FMMM_TILE_OWN 8
每個 thread 用 grid-stride 的方式，最多接 8 個 owner nodes：
for (int i = tid;
i < num_nodes && nown < FMMM_TILE_OWN;
i += nthreads) {
idx[nown] = i;
qx[nown] = pos_x[i];
qy[nown] = pos_y[i];
++nown;
}
這看起來有點反直覺。GPU 不是應該讓 thread 越多越好嗎？
關鍵在 tile 的取得成本。
如果一個 256-thread block 每個 thread 只處理一個 owner，一次 2 KiB 的 tile load 服務 256 個 owners。若每個 thread 擁有 8 個 owners，同一次 tile load 最多服務：
256 threads * 8 owners = 2,048 owners
也就是把 cooperative load 和兩次 barrier 的成本，攤在更多 pair calculations 上。
內層還一次展開四個 partners：
repel(owner, tile[k + 0]);
repel(owner, tile[k + 1]);
repel(owner, tile[k + 2]);
repel(owner, tile[k + 3]);
這讓編譯器更容易交錯獨立指令，隱藏 reciprocal 與 arithmetic latency。
代價是每個 owner都需要自己的 qx、qy、frx、fry、lx、ly。owner 數增加，register pressure 也增加。
如果 compiler 無法把它們都留在 registers，就可能 spill 到 local memory。CUDA 的 local memory 名稱聽起來很近，實際上通常位於 device memory，延遲可能接近 global memory。
所以 OWN = 8 是工程上的 trade-off：
OWN 太小
tile reuse 不足
barrier 成本分攤得不夠
OWN 太大
registers per thread 增加
active warps 可能減少
甚至出現 local-memory spill
最佳值不能只靠直覺。要用實際編譯器 register report 與 Nsight Compute 的 occupancy、local load/store 指標確認。
Naive attraction 很自然地讓一個 thread 處理一條 edge：
thread e = (u, v)
compute force
atomicAdd(force[u], +f)
atomicAdd(force[v], -f)
問題是 high-degree node 會讓很多 edge threads 同時更新同一個 endpoint：
e1
\
e2 ------ hub ------ e3
/
e4
many threads -> atomicAdd(force[hub])
Tile kernel 已經讓一個 thread 擁有 node，因此 attraction 也改成 node-owned traversal。
Graph 先轉成 CSR：
row_ptr[i] .. row_ptr[i + 1]
-> node i 的所有 neighbor IDs
-> 每條 adjacency 的 desired length
owner thread 只累加自己的 attraction：
float f_attr_x = 0.0f;
float f_attr_y = 0.0f;
for (int e = row_ptr[i]; e < row_ptr[i + 1]; ++e) {
int v = nbr_j[e];
float L = nbr_L[e];
// calculate contribution to owner i
f_attr_x += ...;
f_attr_y += ...;
}
因為 accumulator 是 thread-private，這裡不需要 atomicAdd。
代價是 undirected edge 會以兩個 adjacency entries 儲存：
u -> v
v -> u
兩個 endpoints 各自計算自己的 contribution。多一份 adjacency 儲存與讀取，換掉了 atomic contention，也讓結果不再依賴多個 edge threads 的 atomic 執行順序。
這個改動展示一個很實用的 GPU 原則：
改變工作擁有權，
有時比加速 atomic 更有效。
每一輪的所有 forces 必須以同一份 current positions 為輸入：
positions at iteration t
|
v
calculate every displacement
|
grid sync
|
v
update every position
|
grid sync
|
v
positions at iteration t + 1
如果 block 0 算完便先更新座標，而 block 1 還在讀舊座標，同一輪就會混用兩個時間點：
some partners from iteration t
some partners from iteration t + 1
這會把同步 Jacobi-style update 變成與排程順序有關的非同步更新，演算法本身已經改變。
__syncthreads() 只能同步一個 block。Tile kernel 使用 cooperative groups 的：
cg::grid_group grid = cg::this_grid();
grid.sync();
讓整個 grid 在「寫 displacement」與「更新 positions」之間同步。
這也是 cooperative kernel 的限制：所有 blocks 必須能同時 resident。Host 端先查：
cudaOccupancyMaxActiveBlocksPerMultiprocessor(...)
再用：
resident block capacity
= SM count * active blocks per SM
限制 launch 的 block 數量。
這裡的 occupancy 不只是效能數字。它直接決定有多少 cooperative blocks 可以安全參與 grid-wide barrier。
同時，OWN = 8 也有 coverage 限制：
covered nodes <= launched_blocks * 256 threads * 8 owners
因此 backend support check 不能只看 N 與 E 的硬上限；實務上也要測試裝置是否支援 cooperative launch，以及 launch configuration 是否真的覆蓋每個 node。這是把漂亮 kernel 變成可靠 production backend 時很容易漏掉的邊界。
而且這不是假設性的提醒。目前快照中的 tile_force_supported(...) 只檢查：
0 < N <= 200,000
0 <= E <= 2,000,000
如果 cooperative launch 失敗，host code 會改用一個 block 啟動相同 kernel。一個 block 配上 OWN = 8 最多只覆蓋：
1 block * 256 threads * 8 owners = 2,048 nodes
所以 N > 2,048 時，這個 fallback 不能只靠「kernel 有成功 launch」就視為正確。較完整的處理方式應該是：
cooperative launch unavailable
-> select a backend that does not require grid.sync()
-> or split iterations into ordinary kernels at a global barrier
never
-> silently reduce to one block when one block cannot cover N
這也是測試不能只跑主力 GPU 路徑的原因。必須刻意模擬 cooperative launch 不可用，才能驗證 fallback 仍處理全部 nodes。
最直觀的 host-driven loop 會像這樣：
for each iteration:
launch repulsion kernel
launch attraction kernel
launch combine kernel
launch update kernel
synchronize
possibly read convergence data on CPU
Force 計算本身很大時，launch overhead 可能不是主角；但在 multilevel FMMM 的小圖層、很多 seeds 或短 iteration 中，頻繁 launch 與 host synchronization 會變得明顯。
Production path 把 iteration loop 放進 kernel：
for (int it = 1; it <= max_iters; ++it) {
load tiles and compute repulsion;
traverse CSR and compute attraction;
write displacement;
grid.sync();
update positions;
grid.sync();
}
一次呼叫只做：
H2D positions once
H2D topology when needed
one resident multi-iteration kernel
D2H final positions once
其中 positions 仍會在每輪從 global memory 載入 tile；它們不是永遠留在 shared memory。真正省下的是多次 kernel launch、host control，以及每輪返回 Python 的機會。
另一個需要明講的差異是 convergence。Flash kernel 每輪計算平均 displacement，可以在低於 threshold 時提前停止；目前的 tile kernel 雖然從 Python 收到 threshold，device kernel 並沒有使用它，而是固定執行 max_iters：
flash path
stop when average displacement < threshold
current tile path
always run max_iters
固定 iteration count 讓 resident loop 與跨 block 同步更單純，代價是已經收斂的 layout 仍可能多做工作。Benchmark 若只寫「相同 threshold」，兩個 backend 其實沒有相同的停止條件；公平比較應同時報告實際執行輪數。
專案還用可成長的 workspace 保存 device buffers：
capacity enough
-> reuse existing cudaMalloc buffers
capacity insufficient
-> grow capacity and allocate again
這避免在每一輪或每一個相近大小的呼叫中反覆 cudaMalloc / cudaFree。
若 graph 足夠小，全部 positions 可以直接留在一個 block 的 shared memory：
project limits
N <= 768
E <= 4096
Flash kernel 使用一個 block 處理一份 layout，shared state 包括：
s_pos[768]
s_disp[768]
s_last[768]
reduction scratch
它的資料流變成：
global positions
|
| load once
v
shared positions
|
| many force iterations
v
global final positions
這比大型 tiled path 更徹底地減少 position round trips，也只需要 block-level synchronization。
可是單一小圖只啟動一個 block，GPU 上其他 SM 可能全部閒著。要提高整張 GPU 的 throughput，較好的用法是一次處理多份 seeds：
block 0 -> layout seed 0
block 1 -> layout seed 1
block 2 -> layout seed 2
...
這是另一個反直覺結果：
單一 layout 的 latency 已經很低，
GPU utilization 卻可能很差。
增加 batch，
單一 layout 未必更快，
整體 layouts / second 會提高。
因此 production code 需要依 graph size 分 backend：
small graph
-> flash: whole position state in one block
larger exact graph
-> tile: shared position tiles across all resident blocks
very large graph or unsuitable device
-> Barnes-Hut or another fallback
Barnes-Hut 用空間樹近似遠方 node groups，把 exact all-pairs 的工作量往較低方向壓。但它需要建樹與不規則 traversal。對中小型 graph，exact tiled kernel 可能因為資料流規律而比較快；N 很大時，N * N 的算術量才會迫使我們改用近似法。
第一個是 tile 變大。
larger tile
-> more reuse per load
-> more shared memory per block
-> possibly fewer resident blocks
第二個是每個 thread 擁有更多 nodes。
more owners
-> amortize tile load and barriers
-> more registers
-> possibly lower occupancy or spills
第三個是融合更多 iterations。
more fusion
-> fewer launches and host round trips
-> longer kernel lifetime
-> more state stays live
-> less flexibility for convergence checks and fallback
所以調校時不能只問：
shared memory 有沒有用？
應該同時量：
kernel duration
DRAM bytes read/write
shared-memory throughput
registers per thread
local load/store, to detect spills
achieved occupancy
eligible warps per cycle
barrier stall reasons
Shared memory 是用容量與同步換取重用。只看「用了 shared memory」這件事，無法判斷最佳化是否成功。
GPU force kernel 使用 float32，CPU reference 使用 float64。加法順序、reciprocal 指令與平行執行都可能產生小幅差異，因此測試不應要求 bitwise equality。
原專案的 unit tests 使用：
np.testing.assert_allclose(
gpu_force,
cpu_reference,
rtol=1e-4,
atol=1e-5,
)
驗證順序應該分層：
1. formula parity
GPU repulsion / attraction 對 CPU double reference
2. iteration invariants
所有 nodes 有被處理
沒有 NaN / Inf
zero distance 不會除以零
displacement 不超過 max_radius
cooperative fallback 仍覆蓋最後一個 node
3. backend differential test
exact、flash、tile 在小圖上輸出接近
4. LCN admission
quantize to integer grid
repair duplicate / node-on-edge / overlapping-edge violations
official validator passes
5. downstream quality
exact Crossing evaluator 計算 initial K
固定 SA budget 後再比較 final K
這裡不能只檢查 force energy 是否下降。Layout 的 force objective 和 LCN 的 crossing count K 不是同一個目標。
force converged
does not imply
minimum crossing count
專案附有 benchmark_fmmm_cuda.py，會分別輸出：
GPU total time
coarsening time
force compute time
CPU NetworkX spring_layout time
但目前快照沒有保存同一張 GPU 上、同一個 commit、同一組參數的 tiled-vs-naive 原始結果，因此這篇不虛構一個 speedup 數字。
正式比較至少要固定：
GPU model and clock policy
CUDA build flags
graph topology and seed
force formula
iteration count or convergence threshold
warm-up policy
whether H2D / D2H is included
whether coarsening and quantization are included
並同時報兩種範圍：
kernel-only time
回答 tiled force kernel 本身快多少
end-to-end layout time
包含 topology preparation、copies、coarsening、quantization
回答使用者實際等多久
如果拿 NetworkX FR 當 CPU 對照，還要明確寫出兩邊不是完全相同的 force model。這種結果能比較「兩個 layout providers 的耗時與下游品質」，不能歸因成單一 CUDA optimization 的純 speedup。
今天唯一可以不靠 GPU 實測而精確回答的，是特定 mapping 下的 load-request 模型：
tile positions = 256
position bytes = 8
global bytes per block per tile
naive, 1 owner/thread = 512 KiB requested
tiled = 2 KiB requested
theoretical request reduction = 256x
實際 DRAM traffic 與 runtime 必須由 profiler 回答。
Day 23 的核心是演算法增量化：移動一個 node，只重算可能改變的 crossing pairs。
Day 25 的 exact repulsion 無法刪掉 node pairs，因此它展示的是另一條路：
work cannot be removed
-> reduce data movement
-> increase reuse
-> keep accumulators close to compute
-> fuse synchronization boundaries
具體來說：
這也是選 Layout 當 GPU 案例的原因。它很清楚地表現出 GPU 效能不只來自更多 threads；真正的差異常常來自同一筆資料被搬了幾次，以及中間狀態能否留在離運算單元更近的位置。
Day 26 會進入 simulated annealing。Layout 產生一個起點後，SA 會直接針對 LCN objective 提出 move、計算 delta、決定接受或拒絕。第一步先從單條、容易驗證的 Naive SA chain 開始。
iThome鐵人賽
看影片追技術